Three dimensional evolution of differentially rotating magnetized neutron stars 
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We construct a new three-dimensional general relativistic magnetohydrodynamics code, in which 
a fixed mesh refinement technique is implemented. To ensure the divergence- free condition as well 
as the magnetic flux conservation, we employ the method by Balsara [43]. Using this new code, we 
evolve differentially rotating magnetized neutron stars, and find that a magnetically driven outflow 
is launched from the star exhibiting a kink instability. The matter ejection rate and Poynting flux 
are still consistent with our previous finding [3(1 obtained in axisymmetric simulations. 
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I. INTRODUCTION 

Ground-based gravitational- wave detectors, Advanced 
LIGO, Advanced VIRGO, and KAGRA will be in oper- 
ation in the next five years The first observation of 
gravitational waves, thus, will be achieved in the near fu- 
ture. Among their sources of gravitational waves, coales- 
cence of binary neutron stars (BNS) is the most promis- 
ing one, and the detection of gravitational waves from 
them will provide us unique information of strongly grav- 
itational fields and properties of dense nuclear matter. 
The BNS merger is also the potential candidate for the 
progenitor of short- hard gamma-ray bursts 0]. For the 
theoretical studies of the BNS mergers, numerical rela- 
tivity is the unique and robust approach. A number of 
numerical simulations have been performed [Ml since 
the first success in 2000 [2"oT |. 

Magnetic fields could play an important role in BNS 
mergers because the inferred value of the magnetic-field 
strength via the observed spin period P and their time 
derivative P is high as 10 11 — 10 14 G for radio pul- 
sars, of which more than 1800 are known to date [211 ] . 
During the merger process, several mechanisms such as 
compression, magnetic winding, and magnctorotational 
(MRI) instability [22| could amplify their magnetic-field 
strength. This amplified magnetic field could have an 
impact on the dynamics of the mergers because it con- 
tributes to the angular momentum transport and the 
magnetic pressure may modify the structure of the ob- 
jects formed after the merger. Motivated by this expec- 
tation, several groups have implemented the magnetohy- 
drodynamics (MHD) code in the framework of numeri- 
cal relativity [23rl28j . These numerical codes developed 
have been applied to collapse of magnetized hypermassive 
neutron stars (HMNS) [23j2lL magnetized neutron star- 
black hole binary merger j32L l33j ] , evolution of magnetized 
neutron stars 13414361 1 . and magnetorotational collapse of 



massive stellar cores 



In the context of BNS mergers, several groups have as- 
sessed what the role of magnetic fields is during the inspi- 
ral and merger [3814421 ] . Their findings are summarized as 
follows: As long as the magnetic-field strength before the 



merger is not unrcalistically large, e.g., 10 16 -10 17 G, the 
magnetic field does not give a strong impact on the inspi- 
ral dynamics. When the external layers of the two neu- 
tron stars come into contact, the Kelvin-Helmholtz insta- 
bility develops and forms vortexes. Poloidal magnetic- 
field lines are curled by them and generate a toroidal 
field in a short timescale. The saturation point of the 
magnetic-field strength is still under debate because the 
field strength found in Ref. [ZtJ is not as high as found in 
Ref. (43| . If the total mass of BNS is large enough to col- 
lapse to a black hole surrounded by a torus, the mag netic 
field in the torus may be subject to the MRI [41| . On 
the other hand, if the total mass is not large enough, a 
long-lived HMNS [73| is born and the magnetic-field am- 
plification would be realized inside the HMNS. The later 
case has not been explored in detail, because high com- 
putational costs for a longterm well-resolved simulation 
prevent this study. 

The recent measurement of mass for PSR J1614-2230 
(Mji9i 6 _2230 = 1.97 ± O.O4M ) 0] gave a strong con- 
straint on nuclear equations of state (EOS). Together 
with the fact that the canonical observed mass of neutron 
stars is 1.3-1.4Mq, it is natural to infer that a long-lived 
HMNS will be born in the merger of BNS composed of 
neutron stars with the canonical mass Q. This implies 
that it is mandatory to perform a longterm and high- 
resolution simulation of magnetized BNS mergers. 

In Ref. [35[, we have developed a three-dimensional 
general relativistic magnetohydrodynamics (GRMHD) 
code, which has an uni-grid structure with fish-eye co- 
ordinates. The dynamical range of the BNS system is 
quite large spanning from the neutron-star size to the 
wave length of gravitational waves. Thus, we should im- 
plement a mesh refinement technique to save computa- 
tional costs. In the framework of GRMHD, the imple- 
mentation of mesh refinement techniques has been done 
in three methods. For all the approaches, a special care 
for preserving the divergence-free condition of magnetic 
fields is taken. In the first approach, equations for vector 
potentials instead of magnetic fields are solved. In this 
method, any unconstrained interpolations in the refine- 
ment boundaries, where the boundary condition for child 
domains are determined using the data of their parent do- 
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mains, may be allowed for preserving the divergence-free 
condition of magnetic fields [4(J EE El] . In the second ap- 
proach, the hyperbolic divergence-cleaning prescription is 
employed to ensure the divergence- free condition of mag- 
netic fields [H, HU ■ The third one [HI is based on the 
Balsara's constrained-transport scheme, in which a spe- 
cial interpolation scheme in the refinement boundaries 
is mandatory to preserve the divergence-free condition 
and the magnetic-flux conservation [47l E3- We con- 
struct a new GRMHD code employing the third approach 
(modifying the original scheme for the use in the vertex- 
centered grid) to precisely guarantee the divergence-free 
condition and the magnetic-flux conservation. It is wor- 
thy to note that this method is likely to work well also 
in the presence of a black hole. 

As the first application of our new code, we extend 
our previous work in Ref . [36| , in which an axisymmetric 
HMNS with magnetic fields was evolved. In that work, 
we found that a mildly relativistic outflow is driven from 
the HMNS accompanying a strong Poynting flux of mag- 
nitude proportional to B 2 R?Q, (where B, R, and f2 de- 
note the typical magnitudes of the magnetic field, radius, 
and angular velocity of the HMNS) emitted toward the 
direction of the rotational direction. However, it was 
not clear that three-dimensional effects, in particular the 
effect of nonaxisymmetric instabilities such as kink in- 
stability [62j |. would not play a role in this phenomenon. 
For a more physical study, we obviously had to perform 
a three-dimensional simulation. 

The paper is organized as follows: In Sec. [Hi the for- 
mulation to solve Einstein's equations as well as GRMHD 
equations arc briefly summarized. In addition, we briefly 
describe a method to implement the fixed mesh refine- 
ment (FMR) algorithm in particular for magnetic fields, 
and also mention the initial condition and grid setup. In 
Sec. IIIH we present numerical results for the evolution of 
a rapidly rotating magnetized neutron star, focusing on 
the properties of the material and Poynting flux ejected 
from it. Section IIVI is devoted to discussing the impli- 
cation of our numerical results and a summary of this 
paper. In the Appendix, our method for implementing 
the FMR scheme and results for the several standard 
test-bed simulations are shown. Throughout this paper, 
Greek and Latin indices denote the spacetime and spatial 
components, respectively. 



II. FORMULATION, METHOD AND MODEL 

A. Formulation and numerical issue 

We study the evolution of a rapidly rotating magne- 
tized neutron star by a three-dimensional GRMHD sim- 
ulation in the framework of ideal MHD. The formula- 
tion and numerical scheme for solving Einstein's equa- 
tions are the same as in Ref. 



in which one of 
the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) for- 
mulations |49|-|52| is employed, and a fourth-order finite- 



differencing scheme in the spatial direction and a fourth- 
order Runge-Kutta scheme in the time integration are 
implemented. The advection terms in Einstein's evolu- 
tion equations are evaluated with a fourth-order lopsided 
finite-differencing scheme, as, e.g., in Ref. [Hj]. A conser- 
vative shock-capturing scheme is employed to integrate 
GRMHD equations. Specifically, we use a high resolu- 
tion central scheme [54j with the third-order piece-wise 
parabolic interpolation and a steep min-mod limiter. 

We implement a FMR algorithm to our original three- 
dimensional GRMHD code [35[ which has an uni-grid 
structure with fish-eye coordinates [55|. Our FMR 
scheme is essentially the same as an adaptive-mesh re- 
finement (AMR) scheme of SACRA Q, and enables us to 
assign fine grids in the vicinities of neutron stars or black 
holes, while enlarging the computational domain which 
covers a local wave zone of gravitational waves with less 
computational cost. The schemes for solving Einstein's 
equations and hydrodynamics equations (the continuity, 
momentum, and energy equations) are also the same as 
those of SACRA code Q , in which geometric variables and 
fluid variables (density, pressure, internal energy, specific 
momentum, and specific energy) are placed at the vertex- 
centered grids and the grid spacing of a "parent" domain 
is twice as large as that of its "child" domain. Each 
domain is equally composed of (2N + 1, 2N + 1, 2N + 1) 
Cartesian grid zones for [x, y, z), each of which covers the 
interval [— N Axi, NAx{\ for the x-, y- and z-directions 
with Axi being the grid spacing of the Z-th FMR level. 
The label Z varies from 1 (for the coarsest and largest 
domain) to Z max (for the finest and smallest one). The 
prolongation, i.e., interpolation from a "child" domain 
to a "parent" domain, of the geometric and fluid vari- 
ables (i.e., the interpolation of the data of the parent 
domain for determining the data in the child domains) 
in the refinement boundaries are done with a fifth-order 
Lagrange interpolation. Because our grid is located at 
the vertex centers, the restriction procedure, i.e., inter- 
polation from a "parent" domain to a "child" domain, is 
straightforwardly done, by simply copying the data from 
a child domain to its parent domain. 

On the other hand, for integrating the induction equa- 
tion, we need a special care to preserve the divergence- 
free condition of magnetic fields and to guarantee the 
magnetic-flux conservation. For this purp ose, several 
GRMHD codes constructed so far [1^, [13 HM, S3 have im- 
plemented either the constrained-transport (CT) [56[ or 
fiux-CT scheme [57| ■ In the code of implementing FMR 
or AMR algorithm, we are required to take an addition- 
ally special care when performing the prolongation and 
restriction procedures of magnetic fields in the refinement 
boundaries, as described in Refs. AMR- GRMHD 

codes of Refs. [4(1 EH exploited a method of evolving the 
vector potential. This method guarantees the preserva- 
tion of the divergence-free condition for magnetic fields 
avoiding complex interpolation procedures in the refine- 
ment boundaries. However, this method does not guar- 
antee the magnetic-flux conservation in the refinement 
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boundaries, because the magnetic fields are calculated by 
a finite differencing of the vector potential and this pro- 
cedure does not in general guarantee the magnetic-flux 
conservation. 

Alternatively, the AMR-GRMHD code of Ref. [jjl em- 
ploys a hyperbolic divergence-cleaning technique [64j . In 
this scheme, a scalar field is introduced, which is coupled 
to the system of the MHD and induction equations. No 
special prescription is needed for the finite differencing 
when solving GRMHD equations, and a non-zero diver- 
gence of magnetic fields either propagates or damps away 
when they are spuriously excited. However, as mentioned 
in Ref. [45|] j this method is likely to be incompatible with 
the moving puncture method [5lLl52T|. which is commonly 
used to evolve black-hole spacetimes in the BSSN formu- 
lations. 

The AMR-GRMHD code in Ref. 01 implements the 
flux-CT scheme for magnetic fields. In this scheme, both 
the preservation of the divergence-free condition and the 
magnetic-flux conservation are guaranteed in the refine- 
ment boundaries in the machine precision level. How- 
ever, the code of Ref. [24[ is second-order accurate both in 
space and in time. We have developed a new code, which 
is based on the flux-CT scheme (i.e., which can ensure the 
magnetic-flux conservation and the divergence-free con- 
dition in the refinement boundaries), and which is fourth- 
order accurate in time. Following the method described 
in Refs. [47l l48j. we employ a divergence- free-preserving 
interpolation based on WEN05 [58|. In Appendix A, the 
grid structure, the scheme for integrating the induction 
equation, as well as the prolongation/restriction schemes 
of magnetic fields arc described in details. Results for 
several test-bed simulations in Appendix B illustrate the 
reliability of our new code. 



B. Initial condition, density atmosphere and grid 
setup 

Following Refs. [13,|3^], we adopt a rapidly and differ- 
entially rotating neutron star in an axisymmctric equilib- 
rium as the initial conditions. It is a model of the HMNS 
formed after the merger of a BNS. To model the neutron 
star, the following piecewise polytropic EOS composed 
of two pieces is employed: 



P, 



cold 



K 2 p T2 



(p < Pnuc), 
(P > Pnuc)- 



(2-1) 



Here, P and p are the pressure and rest-mass density, re- 
spectively. The specific internal energy, e, is derived as- 
suming the first law of thermodynamics de = —Pd{l/p), 
and this specific internal energy written as a function 
of p is referred to as e co id (i-e., we initially set e = 
£ co id(p))- Following Ref. [30j . the parameters are cho- 
sen to be Ti = 1.3, T 2 = 2.75, K x = 5.16 x 10 14 
cgs, K 2 = K lP l^ T \ and p nuc = 1.8 x 10 14 g/cm 3 . 
This EOS produces spherical neutron stars whose max- 
imum gravitational mass M max (rest mass -Mb. max ) is 



2.OIA/0 (2.32M Q ) and rigidly rotating neutron stars with 
A/max (M b , max ) = 2.27Af Q (2.6OM ). For the rota- 
tional law, we assume the same profile as employed in 
Refs. [10, HI)]. Table U shows the parameters of the dif- 
ferentially rotating neutron-star model which we adopt. 

During the simulations, we use a hybrid EOS as P = 
P co \d + (r t h - l)/o(e - £ C oid) with r t h = IY Our choice of 
r t h may be rather small. We choose this small value to 
focus on the mass ejection from the rotating neutron star 
primarily by the magnetorotational effects suppressing 
shock heating effects. 

A dipole magnetic field is superimposed initially. We 
assume that the axis of the dipole is aligned with the 
rotation axis as in the previous paper (3(| , and write the 
vector potential in the form 



A^ — 



A m 2 



(R 2 + z 2 + ml)*/* 



(2.2) 



where we used the cylindrical coordinate (R,z,ip). zuq 
is set to be 10/3i? e with R e being the equatorial stellar 
radius. Aq determines the magnetic-field strength and 
we adjust this parameter to achieve the maximum field 
strength B to be 4.2 x 10 13 G and 1.7 x 10 14 G. Ac- 
cording to the magnetic-field strength, we refer to these 
models as B13 and B14, respectively. Here, B is de- 
fined by B = \JWb^ where & M is the four-vector of the 
magnetic field in the fluid rest frame. 

We note that there is no reason to believe that the 
dipole axis is aligned with the rotation axis for the HMNS 
formed after a BNS merger. The reason for our choice 
of this simple profile is that the purpose of this paper 
is to compare the results in three-dimensional simula- 
tions with those in the axisymmetric one performed in 
Ref. [36|. If the axes of the dipole and rotation do not 
align with each other, the mechanism for the amplifica- 
tion of the magnetic field and subsequent dynamical pro- 
cess of the system could be significantly modified. We 
will perform more systematic studies varying the axis di- 
rection of the dipole in the future work. 

As discussed in Ref. [3(| , a tenuous density atmosphere 
has to be put outside the neutron star for stably evolv- 
ing magnetically driven outflows. If the atmosphere is 
dense, the outflow and magnetic-field profile are substan- 
tially affected by the inertia of the atmosphere. Thus, we 
have to set the density of the atmosphere to be as low as 
possible. Specifically, we set it as 



Pat 



/at Pn 



(r < 2i? c ), 



/at p max (r/2R c )~ n (r > 2R e ), 



(2.3) 



where p max denotes the maximum rest-mass density of 
the neutron star. We set / a t = 10~ 9 and n = 2. As 
long as we use these values, the magnetic field evolution 
depends only weakly on the atmosphere [36 1. 

Table [TT] summarizes the dipole field strength and grid 
setup. The stars are contained in the numerical domain 
composed of the finest grid resolution. In the typical 
simulations, R c is covered by 80 grid points in the finest 



4 



domain. The side length in one positive direction of the 
finest domain is 1.2i? G . We prepare 8 refinement do- 
mains and in this case, the outer boundary is located at 
« 150i? e - For models B13 and B14, the simulations were 
performed in this grid setup. To confirm the convergence 
of numerical results, we also performed a simulation for 
the low-resolution model B13L, in which R c is covered 
by 60 grid points in the finest domain, while keeping the 
outer boundary at the same position as the high resolu- 
tion model. In these three models, the equatorial plane 
symmetry is imposed. We also performed a simulation 
for model B14F for which no symmetry is assumed. This 
simulation is devoted to exploring the dependence of the 
symmetry effect on the evolution of the magnetized neu- 
tron stars. 



III. NUMERICAL RESULTS 

A. Prediction 

First, we summarize the predicted numerical results 
based on the findings of Ref. (36|, and then, describe 
three-dimensional effects that are not taken into account 
in the previous work [36T |. In the present setup, the mag- 
netic winding due to the presence of differential rotation 
and poloidal magnetic fields will take place and then, a 
strong toroidal field will be generated [z3]. The mag- 
nitude of the toroidal magnetic field increases linearly 
with time during the winding. In particular, a strong 
magnetic field is generated near the rotation axis. After 
the substantial winding, the magnetic pressure associated 
with the strong toroidal field overcomes the gravitational 
binding energy in the vicinity of the neutron-star polar 
surface. Then, a sub- or mildly rclativistic outflow will be 
launched primarily toward the direction perpendicular to 
the equatorial plane. In the outflow, both a matter out- 
flow and a Poynting flux are generated. The magnitudes 
of the luminosity for these would be [36| 



L M ~ 10 48 



L B ~ 10 47 



f Bq \ 2 / R c \ 3 / 



V10 13 Gy V!0 6 cm/ \10 4 rad/s 



10 13 G/ \W 6 cmJ Vl0 4 rad/s 



n 



ergs/s, 
(3.1) 

ergs/s, 
(3.2) 



where Q is the typical magnitude of the angular velocity. 
Here, Lm includes the contribution of the rest-mass en- 
ergy flow, and thus, the luminosity for the kinetic energy 
would be by about two orders of magnitude smaller for 
the outflow velocity ~ 0.1 - 0.2 c. 

The simulations of Ref. (3(| were performed assum- 
ing the axial symmetry. In the nonaxisymmetric case, 
we should consider that the nonaxisymmetric kink insta- 
bility [HH could turn on because the outflow contains a 
strong toroidal field generated by the windin g as a dom- 
inant magnetic-field component. Reference [63j indeed 



showed that the kink instability turns on in a magneti- 
cally driven jet from a black hole-torus system, if it has 
a strong toroidal field. In the following, we will show a 
numerical result which illustrates that the kink instabil- 
ity indeed turns on. The question is how the effect of 
this instability modifies the results of the axisymmetric 
simulations 1361. 



B. Properties of outflow 

Figured] plots the evolution of the electromagnetic en- 
ergy as a function of time for all the models employed in 
this paper. We define the electromagnetic energy by 



En = 



1 

8~~ 



(3.3) 



where 7 is the determinant of the spatial metric and w = 
—n^u^ with n M being a timelike unit vector normal to 
the spatial hypersurface and u M being a four velocity. 
Eb is decomposed into the poloidal component Ep and 
toroidal one Et as 

E ? = Z- I ( bRbR + b z b z )w^/jd 3 x, (3.4) 

87T J 

bn^Wy/^f^x. (3.5) 



E 



1 

8^ 



As expected, the toroidal-field energy for all the models 
increases with time due to the magnetic winding. In a 
relatively short timescale ~ 1 ms (ps 2P c ), the toroidal- 
field energy catches up with the poloidal one, and then, 
it becomes the dominant component. We find the growth 
rate of the toroidal field is consistent with the winding 
mechanism, in which the toroidal field Bt should increase 
in proportional to ~ BpQt with Bp being the poloidal 
magnetic field. 

The toroidal-field energy for model B14F starts de- 
viating from that for model B14 at t ~ 3 ms and the 
growth rate for model B14F is slightly smaller than that 
for model B14. On the other hand, for t > 24 ms, the 
toroidal-field energy for model B14F overcomes that for 
model B14. These facts imply that an asymmetry with 
respect to the equatorial plane comes into the play in the 
amplification process of the magnetic field (see also the 
left panel of Fig. although this effect does not change 
the amplification process qualitatively and significantly. 

For all the models, the poloidal-ficld energy also 
changes with time, and is eventually larger than the ini- 
tial value by a factor of ~ 2 - 10. If the system is ax- 
isymmetric, the poloidal-ficld energy changes only by the 
motion in the meridional plane. Assuming the conserva- 
tion of the magnetic flux and mass, the poloidal field in- 
creases in proportional to p 2 ^, i.e., by compression (see 
e.g., Ref. [37| )■ However, as displayed in the right-panel 
of Fig. [TJ the central density is approximately constant 
during the evolution for all the models. This indicates 
that the increase of the poloidal-ficld energy is not due 
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TABLE I: Physical parameters of a differentially rotating neutron star employed: Gravitational mass 
M, baryon rest mass Mb, central density (maximum density) p max , angular momentum cJ/GM 2 , 
central rotation period P c , and coordinate radius on the equator R e . 

M (M Q ) M b (M e ) p max (g/cm 3 ) cJ/GA/ 2 P c (ins) R c (km) 
2.02 2.23 9.49 x 10 14 0.66 0.48 11.4 



TABLE II: Model parameters and grid setup: Maximum strength for the initial dipole magnetic 
field Bo, the finest grid resolution Aa;z max , the grid point within one refinement domain N, the 
total number of FMR domains imax, the location of the outer boundary Lo along each axis, and the 
assumption for the equatorial plane symmetry. 

Model Bo [G] Ax; max [km] N Z max Lo [km] eq-symmetry 

B14 1.7 x 10 14 0.142 

B13 4.2 x 10 13 0.142 
B13L 4.2 x 10 13 0.190 
B14F 1.7 x 10 14 0.142 



96 8 1740 yes 

96 8 1740 yes 

72 8 1740 yes 

96 8 1740 no 



to the compression, but to a nonaxisymmetric effect. We 
will revisit this point below. We note that the internal 
energy is also approximately constant during the simula- 
tions, which also implies that the density distribution in 
the inner region of the neutron star is approximately sta- 
tionary. Therefore, the magnetic field evolves passively 
with respect to the bulk motion inside the neutron star. 
This is also recognized from the left panel of Fig. [TJ which 
shows that the profile of the curves of Eb for models 
B13 and B14 is quite similar besides the factor deter- 
mined by the ratio of the initial magnetic-field strength. 
By contrast, the magnetic field actively evolves outside 
the neutron star, as described below. This active non- 
linear evolution slightly modifies the scaling relation of 
Eb that might hold between the models of different ini- 
tial magnetic-field strengths. 

Figurcs[2]and|3]plot the snapshots of the magnetic-field 
strength on a meridional plane (x-z plane) at selected 
time slices for models B14 and B14F, respectively. The 
initially dipole field is distorted by the magnetic wind- 
ing and consequently the strong toroidal field is gener- 
ated near the rotation axis (see the central part in the 
middle panels of Figs. [5] and [3]). Then, an outflow is 
driven in particular along the rotation axis (see the mid- 
dle and right panels of Figs. [2] and [3]). The outflow keeps 
blowing for a timescale much longer than the dynamical 

— 1/2 

one ~ Pmax and rotation period P c , because the strong 
toroidal field continues to be generated by the differen- 
tial rotation in the neutron star. The asymmetry with 
respect to the equatorial plane develops for model B14F, 
which causes the less efficient winding as shown in Fig. [1] 
This may be also found by comparing the magnetic-field 
strength in the equatorial plane between two models; for 
model B14, the magnetic-field strength there is weak be- 
cause of the symmetry imposed, whereas it is stronger 
for model B14F. Namely, the winding occurs less coher- 
ently for model B14F. This less coherence is likely to 
stem from the fact that the kink instability turns on in 



a stronger way in the absence of the equatorial plane 
symmetry (see Sec. lIIICj) . Because of this less coherence, 
the toroidal field grows with a longer timescale, and as 
a result, the head speed of the outflow for model B14F 
is slightly slower than that for model B14 (compare the 
right panels of Figs. [5]and|3]). 

Figure [4] plots the profiles of the magnetic field, rest- 
mass density, and z-component of the three velocity along 
the z-axis at several selected time slices for model B14. 
In the vicinity of the stellar polar surface, the strong 
magnetic field is generated and its strength reaches up 
to ~ 10 15 G. This substantially winded magnetic field 
causes a mass stripping if the magnetic pressure over- 
comes the gravitational binding energy density. This is 
approximately written by B 2 /8ir > pGM/H where H is 
the vertical coordinate radius of the neutron star [3fj| . 
This is approximately equivalent to 

V a > V csc , (3.6) 

where v a and v csc denote the Alfven speed and the es- 
cape velocity from the stellar surface. After the substan- 
tial winding, this condition is satisfied for a low-density 
surface region of the neutron star, and mass stripping 
turns on. Figure @] indeed shows that near the stellar 
surface, z = H ~ 1 km, at t = 6.01 ms, the Alfven 
speed ~ \J B 2 /(Airp) is of order the speed of light, c, and 
hence, Eq. (|3.6p is satisfied. After the mass stripping sets 
in, a blast wave is generated, and the shock front propa- 
gates along the z-dircction with its sub-relativistic head 
speed ~ 0.1 - 0.2 c. We infer that differentially rotating 
magnetized neutron stars will drive the outflow as far as 
it is alive. 



C. Kink instability 

A noteworthy new feature that was not able to be 
found in Ref. [36| is that a nonaxisymmetric structure 
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FIG. 2: Snapshots of the magnetic field strength on a meridional plane (x-z plane) at t = ms (left), t = 15.0 ms (center), and 
t = 24.3 ms (right) for model B14. 



of the magnetic-field profile with respect to the rotation 
axis emerges in the outflow; see the middle and right 
panels in Figs. [2] and [3l This implies that a nonaxisym- 
metric instability sets in. This nonaxisymmetric profile is 
found in particular in the vicinity of the rotation axis. As 
already described, the strength of the toroidal field gen- 
erated by the magnetic winding is as large as or larger 



than the poloidal-ficld strength in the region near the 
rotation axis. It is well known that a cylindrical plasma 
column surrounded by toroidal magnetic fields causes the 
kink instability [62} . The situation for the vicinity of the 
rotation axis is quite similar to the cylindrical plasma 
column. This instability is known to turn on if the follow- 
ing Kruskal-Shafranov (KS) instability criterion is satis- 
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FIG. 3: The same as Fig. [2] but at t = 3.11 ms (left), t = 15.0 ms (center), and t = 27.2 ms (right) for model B14F. 
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fied tea: 



> 



2na 
Rq 



(3.7) 



Here, a and Rq are the radius and poloidal extent of 
a cylindrical column, respectively. Because a would be 
smaller than R e and Ro sa z, the condition (|3.7[) can be 
easily satisfied once the toroidal-field strength is compa- 
rable to the poloidal-ficld one. Hence, we infer that the 
magnetic outflow driven from the neutron star is subject 
to the kink instability. 

To determine the dominant mode of the kink instabil- 
ity, we perform a Fourier mode analysis for the toroidal 
field by calculating 



C m (t, R, z ) = J B T (t, R, if, z )e 



°d<f. 



(3.8) 



Here, the spatial hypersurface is sliced for a sequence 
of z = zo(=const) planes on each of which we consider 
rings of radius R and perform the azimuthal integral Q3.8P 
along the rings. Varying the radius of the rings and se- 
lected time, we plot Fig. [5] for zq w 1.9i? c - This figure 
shows that m = 1 mode turns on in particular in the 



vicinity of the rotation axis R < 1 km. We find that 
the ratio |S T |/|B P | « 1 at R « 1 km. The right-hand 
side of the KS condition (|3.7[) is an order of 0.1 with 
a ~ 1 km and Ro ~ 20 km. Therefore, the toroidal- 
field strength comparable to the poloidal-filcd one is large 
enough to induce the kink instability. We find that the 
modes other than the m = 1 mode do not exhibit a re- 
markable growth. This implies that the m = 1 mode 
is the dominant mode, and it does not cause a strong 
non-linear mode coupling. 

Figure [5] plots the Fourier components of the toroidal 
field as functions of time, which arc defined by a volume 
integral as 



D„ 



B T e- imv d 3 x. 



(3.9) 



The left panel of Fig. [6] plots the evolution of D m for 
model B14. This shows that the m = 1 mode is domi- 
nant among the m = 1 - 3 modes as expected from Fig. O 
and its saturation amplitude is at most rj 2 - 3 percents 
of the axisymmetric mode. This again clarifies that the 
non-linearity of the kink instability is not strong enough 
to modify the outflow structure significantly. It should be 
pointed out that this weak non-linearity was also found 
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FIG. 7: Evolution of the matter and electromagnetic energy ejection rates, Lm and Lb, for three models B14, B13, and B14F, 
labeled by "14", "13", and "14F", respectively. 



for the magnetic jet driven from the black holc-torus sys- 
tem in the simulation of Ref . [H} . The features found for 
model B14 also hold for other models: The right panel of 
Fig. [6] shows the evolution of the m = 1 mode for three 
models, showing the weak growth of the m = 1 mode. 



One point to be noted is that the saturation ampli- 
tude for model B14F is slightly larger than that of other 
models. The likely reason is that the absence of the 
equatorial-plane symmetry would enhance the growth of 
the kink instability (in other words, the presence of this 
symmetry would suppress the growth for some channel of 
the kink instability). The stronger excitation of the kink 
instability for model B14F results in a stronger modifi- 
cation of the axisymmetric outflow structure, which we 
found in Fig. [3] 



D. Luminosities 

The magnetic outflow driven from the magnetized neu- 
tron star is accompanied by a large amount of the ejected 
material and electromagnetic waves. We here define the 
luminosities for them by 

L M = -l ded^^—g{T^Y t , (3.10) 

J r— const. 

L B = - I d9dip^(T {cm) ) r t , (3.11) 

J r— const. 

where T^ at ' and T^t are the stress energy tensor for 
the matter and electromagnetic field, respectively, g is 
the determinant of the spacetime metric. We note again 
that Lm includes the contribution of the rest-mass energy 
flow and the contribution only from the kinetic energy is 
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FIG. 8: Evolution of the Lao norm of the divergence-free condition for all the models (left) and the poloidal-field energy _Ep 
and toroidal-field energy Et as functions of time for models B13 and B13L (right). The right panel should be compared with 
the left panel of Fig. Q] 



by about two orders of magnitude smaller than Lm for 
the outflow velocity ~ 0.1 - 0.2 c. 

Figure [7] displays the evolution of these luminosities, 
which are calculated for the extraction radius of r ox w 
420 km. We note that the extraction was performed for 
several radii and we confirmed that their luminosities de- 
pend only weakly on the extraction radii. Comparing 
Figs. 0] and [3 we find that for model B14, the outflow 
front reaches the extraction point at t ~ 12 ms, which 
corresponds to the moment of a quick rise of the lumi- 
nosities. Subsequently, the order of the magnitude of 
the luminosities remains approximately unchanged. For 
other models, the feature of the luminosity curves is es- 
sentially the same. The matter energy flux for model 
B14 attains an order of 10 51 erg/s. Figure [4] shows 
that the outflow density at the extraction point for this 
model is approximately 10 5 g/cm 3 and v z ~ 0.1c. These 
values are consistent with the matter energy flux if we 
assume that the matter is ejected quasispherically, i.e., 
Lm ~ ATtr1^pc 2 v z . For the electromagnetic radiation, 
the luminosities for models B13 and B14 are consistent 
with the scaling relation ([3.21) . which implies that the 
scenario described in Ref. [361 ] is still valid even in the 
presence of the kink instability. On the other hand, the 
electromagnetic luminosity for model B14F is about 10 
times smaller than that for model B14 at the end of the 
simulation. This is because the magnetic winding oc- 
curs less coherently due to a stronger effect of the kink 
instability, as already discussed in Sec. IIII CI However, 
during the longterm evolution, the toroidal-field energy 
for model B14F surpasses that for model B14 (see Fig.Q]) 
and in addition, Fig.[7]shows that the electromagnetic lu- 
minosity for model B14F increases gradually with time. 
This suggests that although the increase timescale of the 
toroidal-field energy and the electromagnetic luminosity 
is rather long, the luminosity for model B14F may be 
eventually as large as that for model B14. 

We find that the scaling relation for the matter 



flux (|3.1[) also holds. We confirm that our result basi- 
cally agrees with that in Ref. (36|. This is because the 
saturation amplitude of the kink instability is not large 
enough to disrupt the coherent toroidal magnetic-field 
profile, as already discussed. 



E. Accuracy check 



Finally, we comment on the reliability of the present 
simulation results. The left panel of Fig.|S]plots the mag- 
nitude of the violation for the divergence-free condition 
and the convergence for the evolution of Eb as functions 
of time. We plot the L^ norm of V • B for models B14, 
B13, and B14F. The divergence-free condition is well sat- 
isfied throughout the simulations, which means that our 
FMR implementation for the magnetic field works well. 

The right panel of Fig. [3] plots the evolution of the 
magnetic-field energy for model B13 with two grid res- 
olutions. The magnetic-field energy for both runs suffi- 
ciently converge until t w 10 ms and the poloidal-field 
energy starts deviating after that, although the behavior 
of the evolution is qualitatively the same. We attribute 
this loss of the convergence to the longterm accumulation 
of the numerical errors such as spurious magnetic recon- 
nections due to the poor resolution. We confirmed that 
the qualitative features of our finding, e.g., the emergence 
of the kink instability and the generation of the outflow, 
are not affected by the grid resolution. Therefore, we 
conclude that the present grid setup is fine enough for 
obtaining scientific results for the evolution of differen- 
tially rotating magnetized neutron stars. 
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IV. DISCUSSION AND SUMMARY 



sorbcd flux at the typical synchrotron frequency is 



A. Discussion 

We here discuss possible electromagnetic signals emit- 
ted by the ejecta from a HMNS formed after the merger 
of BNS, referring to the numerical results in the present 
work. As mentioned in Sec. I, the recent observational 
result of PSR J1614-2230 suggests that the maximum 
mass of spherical neutron stars should be larger than 
1.97 ± O.O4M . This indicates that the EOS of neutron 
stars is stiff, and thus, a long-lived HMNS would be a 
canonical outcome of the BNS merger, if the binaries are 
composed of neutron stars of a canonical mass of 1.3 - 
1.4M© with the total mass ~ 2.7M @. 

Electromagnetic signals should be emitted from the 
ejected material of sub-relativistic motion or ejected elec- 
tromagnetic waves. According to recent studies [65l - [67j 
the ejected material will sweep up the interstellar mat- 
ter and form blast waves. During this process turning 
on the shocked material could generate magnetic 
fields and accelerate particles that emit synchrotron ra- 
diation, for a hypothetical amplification of the electro- 
magnetic field and a hypothetical electron injection. The 
emission will peak when the total swept-up mass ap- 
proaches the ejected mass, because the blast wave begins 
to decelerate according to a Sedov- Taylor's self-similar 
solution. The predicted deceleration time depends on 
the total energy E$ and speed of the ejected material /3oC 
as well as the number density of the interstellar matter 
no for a single velocity outflow as [H| 



2 yrs 



En 



10 49 erg 



1/3 



1 cm 



-1/3 



ft, 
0.2 



-5/3 



(4.1) 

Here, the value of no will depend strongly on the site 
where the merger of BNS happens. If the site is in a 
galactic disk, no would be ~ 1 cm~ 3 , whereas if it is out- 
side a galaxy, the value is much smaller ~ 1CP 3 cm~ 3 . 
By the synchrotron radiation, a radio signal of ~ 0.1 
GHz, which is determined by the self-absorption, could 
be emitted as in the afterglow of gamma-ray bursts [(35[ 
for n ~ 1 cm~ 3 and j3 = 0.2. Then, its luminosity and 
flux would depend on the total kinetic energy. Figure [7] 
indicates that for the model with the initial maximum 
field strength of 10 14 G, the luminosity of the total matter 
energy (including the rest mass, internal, and kinetic en- 
ergies) is ~ 10 50 - 10 51 erg/s. Because the typical speed 
of the ejected material is sub-relativistic with /3q ~ 0.1 - 
0.2, the luminosity of the kinetic energy would be ~ 10 48 
- 10 49 erg/s. The latest numerical-relativity simulations 
indicate that the lifetime of the long-lived HMNS would 
be 0.1 - 1 s [sHHI. Thus, the predicted total kinetic en- 
ergy ejected will be at most ~ 10 49 erg for B = 10 14 G. If 
the magnetic-field strength is smaller than 10 14 G, this 
value is smaller by a factor (Bq/10 14 G) 2 . The unab- 



^2.5 mjy 



E{) 



D 



10 49 erg J V 1 cm- 3 

2 



no V/ 2 fPo_ 
0.2 



300 Mpc 



(4.2) 



and the peak flux at the self- absorption frequency 
is approximately two orders-of-magnitude smaller, i.e., 
O(10)/xJy level. The studies of Ref. (see Table 1 
of it) suggest that the total energy of 10 49 erg for the 
sub-relativistic outflow is not large enough to produce 
a radio signal observable by current and planned radio 
telescopes even for the optimistic value no = 1 cm -3 
(~ 10 50 erg is suggested to be necessary), and our esti- 
mate agrees with their results. Thus, this mass-ejection 
mechanism is unlikely to supply a large amount of the 
mass which generates a sufficiently strong radio signal, 
unless the magnetic-field strength is extremely large, as 
large as that of magnetars for which the field strength 
could be - 10 15 G. 

Alternatively, the authors of Ref. jo6l - [68j (see also 
Ref. [(39[ for the original idea) discuss the signals by the 
radioactive decay of r-process nuclei, which arc produced 
from the neutron-rich material in the outflow, and subse- 
quently decay and emit a signal that may be observable 
by current and planned optical telescopes such as LSST. 
In this scenario, the typical duration of the peak lumi- 
nosity is of order a day or less as [69| 



£pcak 



,0.2 



-1/2 



RL 



1/2 



(4.3) 



and the associated peak luminosity is 



-'peak 



7 x 10 41 erg/s 



/ 



3 x 10- 



Po_ 
0.2 



1/2 



10- 3 M, 



o 



1/2 



(4.4) 



where M* jesc is the total amount of the rest mass ejected 
and / denotes the conversion rate of the energy per rest- 
mass energy in the ejected material through the radioac- 
tive decay process, which is ~ 3 x 10 -6 according to the 
results of |66| . According to Ref. [H, IfJTj , if the total 
ejected mass is > 10~ 3 Af Q , the signal can be detected by 
large optical surveys. Figure [7] indicates that the mass 
ejection rate is ~ 10~ 4 - lO _3 M /s for the maximum 
field strength of 10 14 G. Thus, even if the lifetime of the 
HMNS is 1 s, total amount of the rest mass ejected will 
be ~ 10~ 4 - 1O _3 M for this field strength. Again, un- 
less the magnetic-field strength of the HMNS is extremely 
large (as large as that of magnetars) , the HMNS will not 
eject the material which subsequently can be detected by 
current and planned optical telescopes. 

It should be noted that our estimation is based on 
the magnetically driven outflow. BNS will eject a large 
amount of mass (lO^ 3 - 1O _2 M ) with the velocity 0.2 
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- 0.3c through the dynamical torque that works during 
the merger process [70|. Such a material will be also 
composed primarily of neutrons and produce an amount 
of unstable r-process nuclei, which subsequently decay 
and emit a signal that may be observable by current and 
planned optical telescopes (6(| . The large amount of the 
ejected materials may also contribute to generate radio 
signals interacting with the interstellar matter, as argued 
in Ref. 

Most important finding in the previous [36j and present 
works is that electromagnetic waves are emitted from the 
HMNS together with the mass ejection. This implies that 
even in the absence of the generation of blast waves via 
the interaction with the interstellar material, a strong 
magnetic field is generated. We find that the electro- 
magnetic luminosity is 10 49 - 10 50 erg/s for the model of 
the maximum field strength of 10 14 G. For the hypothet- 
ical lifetime of the HMNS of 0.1 - Is, the total radiated 
energy by electromagnetic waves is ~ 10 49 erg, which is 
larger than the total kinetic energy of the ejected mate- 
rial in our model. Such huge magnetic energy, composed 
primarily of Alfven waves, may be reprocessed efficiently 
to an observable signal as in the solar corona, although 
the mechanism is not clear. To clarify this point, a first- 
principle simulation taking into account detailed physical 
processes, as done, e.g., in Ref. [7l[ for the solar corona 
problem, will be necessary. 

Finally, we should comment on the saturation of mag- 
netic field strength. As shown in Fig. [TJ the magnetic 
field continues to grow at the end of the simulation for 
all the models. The magnetic field strength would satu- 
rate if the magnetic field energy is as large as the kinetic 
energy or the thermal energy of the HMNS. We estimate 
the kinetic energy as ~ 10 53 erg and growth rate of the 
magnetic field energy in Fig. Q] as ~ 10 49 erg/ms, where 
we assume the toroidal magnetic field energy continues 
to grow by the magnetic winding and use model B14F as 
a representative model. Then, the magnetic field would 
saturate at t ~ 100 ms, which is shorter than the lifetime 
of the HMNS ~ 1 s. Therefore, our scaling relation (|3. II) 
(|3.2[) would breakdown after the saturation and observa- 
tional signature might be altered after that. On the other 
hand, if the initial magnetic field of HMNS is not strong, 
e.g., 10 13 G, or the HMNS lifetime is O(0.1)s, the mag- 
netic field strength continues to increase and the HMNS 
would collapse to a BH before the magnetic field saturate 
and Eqs. (f57T ]) -<|3T2 ]l would hold. 



B. Summary 

Using a new GRMHD code implementing the FMR 
algorithm based on the divergence-free interpolation 
scheme of Balsara j47j . we performed numerical simu- 
lations for the evolution of a differentially rotating mag- 
netized neutron star, as an extension of our previous ax- 
isymmetric work [3r| . The magnetic winding mechanism 
generates the strong toroidal field in particular in the 



vicinity of the rotation axis, as in the axisymmetric case. 
Subsequently, Alfven waves propagate primarily toward 
the z-direction along the rotation axis and transport the 
electromagnetic energy. After substantial winding, the 
magnetic pressure overcomes the gravitational binding 
energy density and drives a sub-relativistic outflow. Wc 
found that in this tower-type outflow, the kink insta- 
bility develops due to the presence of a strong toroidal 
magnetic field, and modifies the structure of the outflow 
that would have an axisymmetric structure in the ab- 
sence of this instability. However, this instability satu- 
rates in a relatively small amplitude level, and thus, does 
not significantly modify the profile of the outflow. We 
also confirmed that the scaling relations for the matter 
and Poynting luminosities (|3.ip ~ (|3.2p . originally found in 
Ref. [36l ] . hold in the nonaxisymmetric situation as well. 

As mentioned for several times in this paper, the recent 
observation of PSR J1614-2230 suggests that the maxi- 
mum mass of spherical neutron star has to be larger than 
1.97 ± 0.04M Q , and implies that a long-lived HMNS is 
likely to be a canonical product of the BNS merger if 
the binaries are composed of neutron stars with a canon- 
ical mass of 1.3 - 1.4M @. In the formed HMNS, a 
magnetic field will be amplified not only by the magnetic 
winding but also by the MRI. The MRI could cause an 
efficient angular momentum transport. If the strong self- 
gravity of a HMNS is supported primarily by its rapid 
rotation, the angular momentum transport could induce 
the collapse of the HMNS to a black hole surrounded by 
an accretion torus. (This is not the case if a HMNS is 
supported mainly by a thermal pressure [TT| - ) The black 
hole-torus system formed in this scenario is a promis- 
ing candidate of a central engine of short gamma-ray 
bursts [1^, H3, EH, [llf . We also plan to explore this sce- 
nario in the future. The most self-consistent approach 
for the study of these scenarios is to simulate the merger 
of magnetized BNS taking into account a plausible EOS 
for a long time from the late inspiral to the longterm 
evolution of the formed HMNS. We also plan to perform 
this type of the simulations in the future. 
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Appendix A: Code description 

In this appendix, we describe our integration scheme 
for the induction equation and our implementation of the 
FMR algorithm in details. Following Ref. [HI, we choose 
the three magnetic field B^ = y /jn l/ F* l ' fJ - as the basic 
variable where F* u ^ is the dual of the Faraday tensor. 
This magnetic field is purely spatial in the sense that 
n^B^ = 0. Assuming that the ideal MHD condition 
holds, Maxwell's equation is recasted into the divergence- 
free condition and induction equation: 



d a B a = 

d t B a = d b [(B b v a - B a v% 



(Al) 
(A2) 



where v a = u a /u°. We define the corresponding electric 
field by 



E x = -v y B z 
E y = -v z B x 
E z = —v x B y 



v z B v 
v x B z 
v v B x . 



(A3) 



Here, the electric field is related to the flux of the in- 
duction equation by -B b v a + B a v b = e abc E c with e abc 
being the completely antisymmetric symbol in the flat 
space satisfying e xyz = 1. In the following, we do not 
distinguish E a and E a . 



1. Staggered cell 

The most popular and robust finite volume method 
to integrate ideal MHD equations is the constrained- 
transport (CT) scheme [5(1 . In this scheme, a cell for 
the numerical computation is defined so that fluid vari- 
ables are placed in the cell center. Then, the (surface 
averaged) magnetic field and electric field, which is cal- 
culated from the magnetic field and velocity field, are 
placed on the cell surfaces and cell edges, respectively, to 
preserve the divergence-free condition during the evolu- 
tion of the magnetic field. This prescription is compatible 
with an AMR or FMR implementation of Balsara [47( , if 
the cell-centered grid is employed. For this case, the cell 
surface of a parent domain always agrees with the cell 
surfaces of its child domains, and it becomes straightfor- 
ward to guarantee that the magnetic flux penetrating a 
surface in a parent-domain's cell agrees with the sum of 
the magnetic fluxes penetrating the corresponding sur- 
faces in children-domain's cells. 

However, as mentioned in Sec. Ill Al our code is de- 
signed for the vertex-centered grid, because it is suited 
for integrating Einstein's equations. This implies that 
the magnetic field should be placed at each cell in a dif- 
ferent way so that surfaces of defining the magnetic field 
in a parent domain agree with surfaces of its child do- 
mains. Specifically, we define cells and place the mag- 
netic field, as Fig. IH1 shows. Because the vertex-centered 
grid is employed, the geometric and fluid variables are 



placed at the cell corner (i.e., at the grid). We la- 
bel each grid by {i,j,k) for the child domain and by 
(I, J, K) for the parent domain. Then, the magnetic field 
is placed on the cell surface, labeled by (i±l/2, j±l/2, k), 
± 1/2, fc ± 1/2), and (i ± 1/2, j, k ± 1/2) for the 
child cell. The magnetic field is defined in the range of 
[-(N+l/2)Axi, (N + l/2)Axi\. The CT scheme requires 
the flux (electric field) to be placed at the cell edge, la- 
beled by (i±l/2,j,k), (i,j±l/2,k), and (i,j,k±l/2) for 
the child cell. In the following, big letters such as B x and 
E x denote the components of B a and E a in the parent 
domain, and small letters such as b x and e x do those in 
the child domain. 

The parent cell contains the eight child cells and the 
surfaces of the parent cells always overlap with some cell 
surfaces of the child domain. This grid structure is es- 
sential for guaranteeing the divergence- free condition and 
the magnetic-flux conservation. Table lHll summarizes the 
locations where the magnetic and electric fields are de- 
fined. Then, Eq. (|A2[) is discretized straightforwardly as 



d t {B x ) 



d t (By) I+hJtK+ i 



dt(B z ) I+ i J+ i K 
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(E y )i,j + lK 


Azi 
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{E z ) I+1JK+ i - 
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Axi 




{E x ) I+ i J+1K - 


(E x )i+i.j,k 



Ayi 



(A4) 



where (Axj, Ayi, Azi) denote the grid spacing for (x, y, z) 
in the parent level labeled by I. 

According to Refs. [H, [6lj|, the electric field is com- 
puted at the cell edge by the Lax-Fricdrichs formula, 
given by 



E a 



(E 



x\LL 



(E 



x\LR 



(E 



x\RL 



\RR 



C V I nz 

-2 {Br 



B' L ) 
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B y L ) 



(A5) 



at (7+1/2, J, K). To evaluate the flux, the magnetic field 
defined on the cell surface should be interpolated to the 
cell edge where the electric field is defined (see Fig. [S]). 
This means that the magnetic field has the reconstructed 
right and left state. According to the prescription of 
the central scheme [E3| . we need an offset based on the 
characteristic speed in the flux. In the above equation, 
{E X ) RL represents the reconstructed right state in the 
y-direction and left state in the z-dircction. The other 
symbols (E X ) LL , {E X ) LR , and (E X ) RR are interpreted 
in the similar way. (B V R ,B V L ) and (B Z R ,B Z L ) also denote 
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the right and left state of B y and B z reconstructed. c y 
and c z are the characteristic speeds in the prescription 
of the upwind flux construction and calculated at cell 
edges with the interpolated variables. We simply calcu- 
late these quantities by averaging: 



( C y)i+±,J,K 

(,C z ) I+ l JK 



(v V )l,J,K + (v V )l+l,J,K 

2 

(v Z )l,J,K + {v Z )l+l,J : K 



(A6) 



The formula for E y (E z ) is obtained from Eq. (|A"5)l by 
the permutation of the indices x — > y, y — > z, and z — > x 
(x z, y —> x, and z —t y). 

For solving GRMHD equations and Einstein's equa- 
tions, one needs the magnetic field defined at {I,J,K). 
This is done by a simple averaging as 



(B x ) 



(B y : 



I.J.K 



I.J.K 



(B z )i, 



J,K =-. 



( B ' T )/,./+i.A'-i + (B x ) I:J _i iK _i 

(B y ) I+hJtK+i + {By) I+hJ>K _i 
+ {B y ) I ^ K+h + {By) I _ i ^ K _, 

1 (B z )i+i,j+ij + (B z )j_i J+ i K 
+ (B z ) I+ i ! ji tK + (B z ) I _i ! j_ hK 



2. FMR implementation 



For the assignment method of the variables in the 
interior of a cell described in Sec. IA 11 we exploit the 
divergence-free reconstruction scheme in the refinement 
boundaries of the FMR algorithm, following Refs. [13, 
|48| . First, we review how to reconstruct the magnetic 
field in the whole region of a cell in this scheme. 

Consider a cell defined for a domain composed of x £ 
[xi,x I+ i], y G [yj,yj+i], and z £ [z Kl z K+1 ] at a FMR 
level I. In the first step, the magnetic-field profile on 
the cell surfaces are reconstructed. When we design a 
third-order accurate code for the spatial direction, the 
profile of the magnetic field, say B x , on the surface at 
x = xi,y £ [yj,y J+ i],z £ [z K ,z K+1 ] should be written 
as 

B x {x = x I ,y,z) = B% + B x Pi{y) + B x Pi(z) 

+ B* y P 2 {y) + B^P 1 (y)P 1 {z) + B* z P 2 (z), (A8) 

where Pi(y) = y - y J+1 / 2 and P 2 (y) = 
{y - y.J+1/2) 2 - Ay, 2 /12. We employ the WEN05 
scheme to obtain the coefficients Bq, By, B x , B x y , B x z , 
and B x z [48[. In this scheme, we first consider the one- 
dimensional reconstruction problem in a zone centered at 
y = y ,7+1/2, taking into account five neighboring variables 



{(B x )j- 3/2 , (B x )j- }/2 , (B x )j+i/2, (B x )j +3/2 , {B x ) J+5/2 }, 
where we omit the index / and K. Then, a third-order 
reconstruction over the zone centered at y.j+i/ 2 can be 
carried out by using three stencils Si, S 2 and £3 that 
rely on the variables {(B x )j_ 3 / 2 , (B x )j_ 1/2 , {B x ) J+ i /2 }, 
{(B x )j-i/2, (B x ) J+ i /2 , (B x ) J+3/2 }, and 
{(B x ) J+ i /2 , {B x ) J+3/2 , {B x ) J+5/2 }, respectively. Be- 
cause the reconstructed polynomial in the y-direction 
has the form 

B x (y) = B x + B x Pi(y) + B x yy P 2 (y), 

we should calculate B x and B x y for the each stencil in 
the following manner; 

/ dx\(i) 3(B x )j+i/ 2 - 4{B x )j_i /2 + (B x )j_ 3/2 

W = 2AS 

fR1 ^(i) (B x ) J+ i /2 - 2{B x )j_i /2 + {B x )j_ 3/2 
{Byy) ~ 2A^ (A9) 

for the stencil Si, 

{B x ) J+3 / 2 - (B x )j_i /2 



(B X )M 



(B 



v( 2 ) _ (B x )j+ 3 /2 - 2(B x ) J+ i/ 2 + {B x )j_i/ 

> ~ OA.,2 



2Ay? 



(A10) 



(A7) for the stencil S 2 , and 



mx \(3) -(B x )j+5/2 + 4{B x ) J+3/2 - 3(B x ) J+ i /2 
W = 2A^ 



(B 



, x n(3) = (B x )j+5/2 - 2(B x ) J+3 / 2 + (B x ) J+ i/ 2 
yy) ~ 2 Ay? 



(All) 



for the stencil S 3 , respectively. According to the pre- 
scription in Ref. [58|, we calculate the weig ht cjW for 
each stencil with k — 1, 2, and 3. Then, we evaluate the 
coefficients as 

B* = u< 1 \B*)W + ^(B X )^ + ^(B X )^ 

Byy = J'HB^M +JV(B X y )V +JV(B X yy )^. 

(A12) 

The weight u/ fe ) is reduced to be nearly zero if the stencil 
k contains a discontinuity, while, for the smooth profile, 
it is reduced to be the optimal weight, with which the 
right-hand side of Eq. (|A12[) can be a fifth-order accurate 
expression of the derivative. The coefficients B x and B x z 

as well as the cross term B x „ are obtained in the similar 

y ~ 

way. 

Essentially the same procedure is applied to the surface 
at x — xi+i. To reconstruct B y (B z ), the permutation 
rule of x y, y — > z, and z — ¥ x (x — > z, y — > x, and 
z — > y) should be simply applied to. 
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Second, wo reconstruct the magnetic field in the inte- 
rior of the cell, for which the third-order accurate form 
is 

Bf (x, y,z) = a x + a x x P 1 (x) + a^P, (y) + a x z P, (z) 

+ a x xx P 2 (x) + a% y P 1 (x)P 1 (y) + a^P^P^z) 

+ a x yy P 2 (y) + a? y J\{y)P x {z) + a x z P 2 (z) 

+ <xM x ) + <'',,„ I 'r'-) Pin: + a x xxz P 2 {x)P x {x 

+ a x xyy P 1 (x)P 2 (y) + a x zz Px(x)P 2 (z) 

+ a x xyz P 1 (x)P 1 (y)P 1 (z), (A13) 

where P 3 (x) = (x-x I+1 / 2 ) 3 -3(x-x I+1 / 2 )Axf /20. The 

expression for B y (x,y,z) (B (x,y, z)) is also obtained 
from the permutation of x — > y, y — > z, and z — > x 
(x — > z, y — > x, and z — > y). Hence, we have to de- 
termine in total 48 unknown coefficients. Imposing that 
d a B a = holds everywhere inside the cell, we obtain 10 
algebraic equations. Furthermore, we require that the 
profile for the interior of the cell matches that on the 
cell surfaces. Then, 36 algebraic equations are obtained. 
These 46 algebraic equations are not independent, i.e., 
one of them can be derived from the others. This implies 
that there are three degrees of freedom. We fix these de- 
grees of freedom in the following manner: Consider one 
of the algebraic equations 



Finally, using the algebraic equation (|A13[) that holds 
in the whole interior of the parent cell, the magnetic 
fields in the eight cells of the child domain I + 1, con- 
tained within the parent cell Z, are reconstructed. Be- 
cause the algebraic form of the magnetic field satisfies 
the divergence-free condition, the magnetic field in the 
child cells thus determined satisfies this condition auto- 
matically. 

The restriction of the magnetic field from the child cells 
to their parent cell is done at specific time step levels: 
We choose the time-step levels in the FMR algorithm 
following Refs. @,Hl|. Specifically, the time step intervals 
for each FMR level is chosen by 



At, 



c C flAx Zc for 1 < I < l c , 
ccflAxi for l c < I < 

' ma 



(A15) 



where ccfl is the Courant number m 0.4 - 0.5. Namely, 
for the coarser levels with I < l c , the time step intervals 
are chosen to be identical while it is chosen to be propor- 
tional to the grid spacing for I > l c (see Fig.[TU]). In this 
setup, the restriction is done when the time slice of the 
child domain agrees with that of the corresponding par- 
ent domain (see the time step level n + 1 in Fig. [TU|) . The 
simplest form for the restriction would be (cf. Fig. [9]) 



a 4- a y 

^xxy 1 xyy 



(A14) 



where a v xyy 



( a xyz) 1S a coefficient in 
B y (x,y,z) (B z (x,y,z)) in the analogy with Eq. (|A13|) . 
a xyz can be determined by the matching at the cell 
surface and in the interior, and there is no equation 
to determine a* v and a xyy other than Eq. (|A14|) . We 
follow Ref. [47| to determine these coefficients. By 
minimizing the magnetic energy involved in the cell with 
respect to a xxy and a% yy , we obtain 



xxy 



xyy 



(see Ref. [47| in details). The same procedure is applied 
to the coefficients with the permutation of x — > y, y — > z, 
and z — > x, and x — > z, y — > x, and z — > y. As a result, 
three degrees of freedom are fixed. 



(B x ) I>J+ i >K+ i - -(V J+ i fc+ i +b x i j+ 3 k+ i 



(A16) 



However, this cannot be employed, because the 
divergence-free condition is not satisfied in the parent 
level: Note that the divergence-free condition for B a is 
preserved if the flux calculated from the electric field E a 
is used to integrate the induction equation. However, in 
the restriction (|A16[) . e v and e z , instead of E v and E z , 
are used to update B x (cf. Fig. [9]). For such cases, simple 
restriction schemes in general do not work well. 

Thus, following Ref. [47[, we add a correction in addi- 
tion to the "zeroth-order" restriction (|A16[) . to preserve 
the divergence- free condition of the magnetic field. Ref- 
erence 
I > L 



47J proposed to use the following restriction: For 
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■Am) 



Am) 



(B«) I+hJ!lcH H- (B*) 7+ i iJ)Jf+ i + £ Ai| m) - 

m— 1 

.i + £ At (m 



Ax; 



Ax z _ 



4 '+1 



>(m) 



/+i,J+l,A-- 



(•B 2 )/+i,j+i a- -> (B z ) /+ i J+ i x 



m— 1 



(m) 



Ax z 



Aa; i+ i 



l—< 1 At, 4 Z- -< '+ 1 



>(m) 



Az* 4 

m— 1 

H^u i 

+ 4 



Axi+i 



j+4 -k+2 



Ax ;+1 



(A17) 



where the time step intervals of the Runge-Kutta inte- 
gration are defined as follows: Atj, = AtJ, 4 ^ = Afy'/6, 
A4 2) = At; ( , 3) = A^//3 with i' = I, l+l, AtlJj = Afg?! = 
At, +1 /6, and Atg.\ = Atffj = Af, +1 /3. (E»)<™>, 
(E z Y m \ (e y Y m \ and (e z )( m ) denote the electric-field 
components at sub-step levels, m, of the Runge-Kutta 
integration (see Fig. HU)) . For I < l c for which the time 
step intervals are identical, the third term of the right- 
hand sides of (|A17[) should be replaced from ^ m=1 to 
J2m=i- The similar procedure is applied for B y and B z 
by the permutation of the indices. These prescriptions 
guarantee both the magnetic-flux conservation and the 
preservation of the divergence-free condition. 



Appendix B: Code tests 
1. One-dimensional tests 

We here report the results for one-dimensional MHD 
tests in the Minkowski spacetime, proposed in Ref. [59| . 
The initial data of the tests are summarized in Table IIVI 
For all the initial data, a discontinuity is present at x = 0, 
and the left (x < 0) and right (x > 0) states are com- 
posed of uniform profiles. For all the cases, the T-law 
EOS with r = 4/3 is adopted. We performed these tests 
in a three-dimensional code assuming that all the quanti- 
ties are uniformly distributed in the y- and z-directions. 
The divergence-free condition of the magnetic field im- 
plies that B x is a constant along the x direction. Be- 
cause B x does not evolve in this test, the divergence- free 
condition is automatically satisfied. 

The purpose of this section is to demonstrate that our 
FMR code works well even in the presence of discon- 
tinuities and shock waves across the refinement bound- 
aries. To do this, we prepare a computational region com- 
posed of two FMR domains, for which the grid point and 
resolution are summarized in Table IIVI The simulations 



were terminated at <fi n when a discontinuity or waves go 
through the refinement boundary. Figures [IT] - [13] plot 
the snapshots of p and B y at t = ifi n . Numerical solu- 
tions in the FMR domains 1 and 2 are plotted together 
with the red-plus and green-circle symbols. The solutions 
in both domains agree approximately with the analytic 
solutions, except for a spurious small modulation for p in 
the region 1.0 < x < 1.7 in the slow shock test and a small 
bump around x = — 1 in the switch-off test. As reported 
in Ref. [28[, these errors are generated by the initially 
discontinuous data at x — irrespective of grid reso- 
lutions and presence of the FMR refinement boundary. 
The amplitudes of those errors are known to gradually 
decrease with improving the grid resolution. The point 
to be stressed is that the discontinuities and waves pass 
thorough from the domain 2 to domain 1 without any 
problems. These results validate the implementation of 
staggered magnetic fields, prolongation, and restriction 
described in the previous section. 

2. Two-dimensional tests 

We also performed a two-dimensional cylindrically ro- 
tating disk test proposed in Ref. All the variables 
are assumed to be functions only of x and y in this test, 
although the the simulation was performed by a three- 
dimensional code. The employed initial condition was 

(p, P, B x /V^, B y /V^, v x , v v ) 

= ( (10, 1, 1, 0, -uy, lox) for (y ^TF < 0.1) 
[(1,1,1,0,0,0) for (^x 2 +y 2 > 0.1) 

with u) — 9.95. In this test, we adopted the T = 5/3 
EOS and prepared three FMR domains which are com- 
posed of squares with the regions 6 [—0.512,0.512] 
for domain 1, x, y £ [—0.256,0.256] for domain 2, and 
x, y E [—0.128, 0.128] for domain 3. The finest grid spac- 
ing and grid points were chosen to be 0.001 and 128, re- 
spectively, which is equivalent to the middle resolution in 
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TABLE III: Grid points where the geometrical and fluid variables, the magnetic field, and the 
electric field are defined, respectively. 

Metric and fluid variables (/, J, K) 

B x (7,J+i,^+i) 
B y (I+±,J,K + ±) 
B~ {I+\,J+\,K) 
E x (I+±,J,K) 
E y {I,J+\,K) 
{I,J,K+ 1 -) 



Child domain (i + 1 ) Parent domain 0) 




FIG. 9: Schematic picture for the structure of the cells in our FMR algorithm together with the assigned locations for the 
magnetic field in a child domain (left) and in a parent domain (right). Geometrical and fluid variables are defined at the cell 
corner, magnetic field on the cell surface, and flux at the cell edge, respectively. 



Ref. [45| . Comparing the result in Ref. [45| , we found that the bottom-right panel of Fig. IT4| . Therefore, we confirm 
all the quantities were well reproduced. The divergence- that our FMR implementation works well, 
free condition was also satisfied with a high precision (see 
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FIG. 10: Schematic picture of the time integration scheme for a child domain of label I + 1 (left) and parent domain of label I 

from n to n + 1 time slice for the domain I > l c . Attached small arrows with numbers denote the sub-step for the fourth-order 
Runge-Kutta integration. 



TABLE IV: Initial data and grid setup for ID MHD tests. Initial states are separated in the left 
(x < 0) and right (a; > 0) state with a discontinuity at x = 0. 2N + 1 and Z max denote the number 
of the grid point covering the interval [— NAxi, NAxi] in a FMR domain and the number of FMR 
domains, respectively. 



Test 


Left state (a; < 0) 


Right state (x > 0) 




Aa; !max 


N 


^max 


Fast Shock 
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B7^/4^ = (10.0, 18.28, 0.0) B a /V4w = (10.0, 14.49, 0.0) 
P = 10.0, p = 1.0 P = 55.36, p = 3.323 


0.005 


204 


2 


Switch-off 
Fast Rarefaction 
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1.8 
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204 
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Switch-on 
Fast Rarefaction 
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B a /^ = (1.0,1.022,0.0) 
P = 0.1, p= 1.78 x 10" 3 
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204 
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FIG. 11: Snapshots of rest-mass density and y-component of the magnetic field in ID fast shock (left) and slow shock (right) 
problems. Numerical solutions are plotted both for the FMR domain 1 (red-cross symbol) and domain 2 (green-circle symbol). 
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FIG. 12: The same as Fig. [TT] but for the ID switch-off and switch-on test problems. 




FIG. 13: The same as Fig. [TT] but for the shock tube 2 tests. For the ID shock tube 1, we plot the four velocity weighted by 
the enthalpy instead of the magnetic field because B y = in this simulation. 
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FIG. 14: The results ol a cylindrically rotating disk problem: Rest-mass density p (top-left), pressure (top-right), magnetic 
field strength b 2 (bottom-left), and divergence- free condition (bottom-right) on the x-y plane at t — 0.4. 
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